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Molecular interactions are wired in a fascinating way resulting in complex behavior of biological 
systems. Theoretical modeling provides a useful framework for understanding the dynamics and the 
function of such networks. The complexity of the biological networks calls for conceptual tools that 
manage the combinatorial explosion of the set of possible interactions. A suitable conceptual tool to 
attack complexity is compositionality, already successfully used in the process algebra field to model 
computer systems. We rely on the BlenX programming language, originated by the beta-binders pro- 
cess calculus, to specify and simulate high-level descriptions of biological circuits. The Gillespie's 
stochastic framework of BlenX requires the decomposition of phenomenological functions into basic 
elementary reactions. Systematic unpacking of complex reaction mechanisms into BlenX templates 
is shown in this study. The estimation/derivation of missing parameters and the challenges emerg- 
ing from compositional model building in stochastic process algebras are discussed. A biological 
example on circadian clock is presented as a case study of BlenX compositionality. 



1 Introduction 



Computational systems biology is a novel approach to understand how biological systems are orches- 
trated altogether LL(2J. Biology, physics, computer science, systems theory and mathematics have 
joined to propel a multidisciplinary research that provides tools for the analysis of biological stud- 
ies. Particularly, stochastic approaches are becoming more and more popular as novel experimental 
techniques - such as quantitative flow cytometry |3| and fluorescence microscopy [4J - provide single 
level measurements of cell physiology. While the average behavior of a cell population has been de- 
scribed by continuous modeling approaches (e.g. with Ordinary Differential Equations, ODEs) [5 1 from 
a long time, single cells are analyzed in a stochastic framework as fluctuations may have a significant 
effect on the physiology of the cell. The influence of noise also in gene expression and signal trans- 
duction processes have been shown to be important by both theoretical and experimental approaches 
|[6ll3|8l0|TOllIIl|T2l|T3|l4^ Computer science through the discipline of Algorithmic Systems Biology 
|[T5l enables this research. Here we concentrate on a class of formal languages, namely stochastic process 
algebras that have been used for interacting and distributed systems and serve as a modeling framework 
for biological systems as well (for a survey see |16|). Starting from stochastic TT-calculus |[I71[T3| and 
passing through beta-binders [[T9l, a real progranmiing language has been designed specifically to model 
and simulate biological systems: BlenX [201. 

BlenX allows the modeler to create boxes that represent the interacting biological entities (proteins, 
genes, etc.) (Figure [T]). Boxes have an internal program (or internal behavior) describing their behavior 
and a set of typed interfaces describing their interaction capabilities. The interaction sites on boxes are 
called binders. The set of binders {interface) and the internal program of the box drive its behavior. 
The behavior of the biological system is given by the ordered sequence of actions and reactions that the 
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BlenX-based compositional modeling of complex reaction mechanisms 



a: A b: B . . . z:Z 



Internal Program 
[P] 



BOX 



let BOX : bproc = #(a:0,A), #(b:0,B),..., #(z:0,Z) [ P.nil ] ; 



Figure 1: Abstraction of biological entities in Blenx. 



program can perform leading to the biochemical interactions between the elements. Actions for instance 
can occur when binders "sense" signals (receive an input) and propagate signals (send an output) and 
the internal structure codifies for the mechanism that transforms an input signal into a change of the 
box (e.g. activation (unhide or expose), deactivation (hide) or changing the type of a binder). Events 
specify statements to be executed with a rate and/or when some conditions are satisfied. Boxes are 
able to born (when a new box is synthesized) and to die (when a box is deleted) as biological entities 
(Table [T]). Boxes are merged or split upon different conditions (join and split events). They can also 
form complexes through their binders and dissociate depending on the state of the overall system. For 
a more comprehensive description of BlenX we would refer the reader to EOl . Similar process calculi 
initiatives have been also developed for building biological models and performing stochastic simulations 
(i.e. stochastic TT-calculus HHIISJ, BioAmbients [21J, Brane Calculi [22], CCS-R |23|, /c-calculus L24il . 
Bio-PEPA |25|. 



Table 1: Basic primitives of BlenX. A declaration of a binder or a box can be nameless or named, e.g. 
it can have an Id or boxld, respectively. Events are used to express actions that are enabled by global 
conditions, expressed by cond. Actions are that processes can perform. 



EVENTS: synthesis of a box 


when(cond or rate) new(Decimal) 


degradation of a box 


when(cond or rate) delete(Decimal) 


division of a box 


when(cond or rate) split(boxId,boxId) 


join of a box 


when(cond or rate) }oin(boxId) 


ACTIONS: exchange an input/output signal (communication) 


a?().process / b!().process 


expose a binder 


expose(rate, Id: rate, Id) 


hide / show the binder of a box 


hide(ratejd) 1 unhide(ra^^,/J) 


change the type of a binder 


c\v(rate,Id,Id) 



Composition is the first challenge that the modelers should deal with. A key innovative aspect of BlenX 
is the ability to model the reactions between components simply by listing their affinity and without the 
need of programming all the possible interactions. The BlenX framework allows the user to build systems 
by fixing each reaction of the network (also called as bottom-up approach) or gives opportunity to handle 
abstractions as well (such as a top-down approach). After specifying the system, the BlenX program is 
executed with the Gillespie SSA algorithm ll26l . The reactions occurring in the system are defined by rate 
dependent functions that are crucial for the reaction propensities of the stochastic model. Rate functions 
are associated to actions and events of boxes, and those rates can be determined by the mass action 
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kinetic law. A crucial point of biological models built upon mathematical formalisms is the additional 
presence of the complex mathematical functions (e.g. Michaelis-Menten kinetics |27|, Hill function |28|, 
etc.) that have been empirically developed through several assumptions in order to reduce the size of the 
system. These abstractions simplify the system leading to a decrease in the required computational power 
for calculation. Furthermore, modelers often turn to these phenomenological functions to describe the 
observed behavior of a system without knowing all its details, such as multi-step reactions are often 
assumed to happen at the same time in cooperative reaction schemes [29|. Complex rate functions raise 
several problems in stochastic process algebra approaches. 

First, these frameworks are based on Gillespie's stochastic algorithm which considers only elemen- 
tary reactions, while biological models often deal with nonlinear terms in the deterministic framework. 
Nonlinearity is known to serve oscillations in several periodic biological systems 1 30 J or multsistability 
in others lISTIl , giving an important role for these mathematical formulas in simple models. 

Secondly, when modeling has to deal with several assumptions of phenomenological modules, the 
freedom of compositionality is reduced. The hidden parts of the modules may contain important linkage 
between models that are chosen to be merged together. Several biological models are built each day 
making compositionality to be one of the most important key features of process algebra that offers 
systematic modeling of complex biological networks. Compositionality has been addressed as an issue 
of model-construction from elementary reactions with the basic operators |[32l|33l, as the translation of 
one approach to another f34\, or as the combination of different types of models (ODEs with process 
algebras. Boolean, hybrid models) |[35l , but the compositionality of complex rate functions has been 
attracted less attention. 

When non-elementary reactions occur and compound mathematical formulas are used, the direct 
translation of mathematical terms into the stochastic context is a well-liked approach. Usage of these 
general functions for calculating the rate of a reaction is also possible in BlenX ||20l|36l. However, these 
implementations have been pointed out by several authors to be incompatible for some cases |[25l[37l[38l 
[39l|40l|43l, thus modelers have to be careful with them. Stochastic modeling of complex functions is 
only an approximation and assumptions have to be handled globally. Thus the BlenX framework calls for 
a semi-automatic method of describing these complex rate functions with intermediate steps (we refer 
as an "unpacking" mechanism) not only owing to ease the compositional programming process, but to 
provide a correct (and generalized) way of stochastic simulations. 

"Unpacked" version of these modules, representing complex reactions, has to be available as an 
option for merging them properly when a hidden link is becoming to play a role in the whole system. 
Usage of decomposed modules containing only elementary steps lengthens simulation time, but provides 
formal grounds for composition. After merging the modules properly, reduction and abstraction tech- 
niques can be applied to the overall system. It is our belief that "unpacked" modules could be a natural 
way of modeling biological processes by connecting small models into large ones. Modern modeling 
techniques should provide a framework where the model-building process can be carried out correctly 
(assumptions are taken into account), while it is straightforward and is also suitable for the stochastic 
simulations. The main goal of this article is to present a possible extension of process algebra tools to 
carry out compositional modeling in a proper and an easy way. 

In this study we concentrate on BlenX. We propose a novel framework for constructing models by 
previously coded model fragments stored in a template-library. Building large models starting from the 
basic principles of a process calculus language is time consuming and error-prone, thus a template library 
should ease the work of programming as it happens in all the fields affected by computer technology. 
Important motifs with special dynamics have been already proposed in biology ||44l|45l and in computer 
science ll46l , but they rely on compound phenomenological description and mathematical assumptions. 
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By using BlenX primitives, subsystems (called motifs) are defined, templates are built and decomposition 
of complex formulas into elementary steps is achieved. Furthermore, compositional modeling is applied 
to an oscillating biological system (called circadian clock L47J ). We show that template-based modeling 
in BlenX improves compositionality. 

The paper is structured as follows. Section 2 and 3 show the stochastic decomposition of Michaelis- 
Menten kinetics and of Hill functions, respectively. Section 4 uses the stochastic decompositions to 
compositionally build a model of the circadian clock. We conclude the paper with a discussion of the 
results and of the future research directions. 

2 Michaelis-Menten module 

Most of the biochemical reactions require catalytic molecules (called enzymes) that increase the rate of 
the reactions. In enzymatic reactions, the molecules at the beginning of the process are called substrates 
(5), and the enzyme selectively converts them into products (P). The kinetic description of such systems 
was expressed by L. Michaelis and M.L. Menten [27]. The derived equation of their results by (referred 
as Michaelis-Menten kinetics) are widely used in kinetic modeling. The scheme of a one-substrate- 
one-product reaction (with one active site) is presented in Table [2] The catalytic step is supposed to be 
irreversible and the rates of the reactions are given by the law of mass action. 



Table 2: Steps of the enzymatic reaction. 







(1) 


Rate of ES complex formation 


ki E S 


ki 


-ES^E + P 








E + S^= 


(2) 


Dissociation rate of ES complex 


k2-ES 


ki 












(3) 


Production rate of P 


k3-ES 



As enzymes are specific to their substrates and the formation of the enzyme- substrate complex (ES) is 
assumed to be relatively fast, the equilibrium is reached rapidly and the production of P is the limiting 
step in the overall system. Therefore the ES complex is assumed to be stable, thus the change of its 
concentration approaches zero (quasi- steady- state assumption (QSSA)). Another important assumption 
is that the concentration of the substrate highly exceeds the one of the enzyme ([5] » [Etot]). When 
a critical substrate concentration is reached, the enzyme is saturated and additional amount of substrate 
does not influence the velocity of the reaction; it is already maximal (Vmax)- If the last reaction is assumed 
to be irreversible and all the previously mentioned statements are valid, the rate of the substrate turnover 
to product is approximated by 

V Vmax • K^]^ whove Vmax h ' [Efot], (^3 + h)/^! and [Efot] ^E^ES. 
This equation provides a complex rate function assuming a single step reaction: E -\-S ^ E + P 

The Michaelis-Menten rate law is often found to be a good approximation to describe enzymatic reac- 
tions. The rates of complex formation and its dissociation are rarely available in biological systems, 
while the key parameters {Vmax and Km) of a Michaelis-Menten reaction might be easily determined from 
measured data through linear graphical representations (e.g. LineweaverBurk plot, HanesWoolf plot, 
EadieHofstee diagram) or by nonlinear regression methods. 
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In the next subsection, we provide a brief description of how to code enzymatic reations in BlenX 
with elementary steps and we give a hint how to search for unknown parameters in a Michaelis-Menten 
module. 



2.1 Decomposition of the module 

The use of Gillespie's stochastic algorithm requires elementary steps instead of complex rate functions 
in a model. Decomposition of the Michaelis-Menten rate law into elementary reactions may lead to 
crucial changes in the system's behavior as nonlinearity may disappear if assumptions are inconsistent 
about enzyme- substrate complexes ||48l|49l. Compositional model building should carefully handle the 
enzyme molecules hidden in the QSSA. 



x:S 



communication 

-^ rV*.; |----| (Ki) \""\ 



x:P 



x!().ch(s,P).nil 
e: E 

J 

E 

rep e?().nil 



reverse binding'' 



^ ' ctiange ^ 



P 

nil 



LIU J 



(infinite) decomposition 



Figure 2: BlenX representation of the Michaelis-Menten kinetics. 



The Michaelis-Menten module can be implemented easily into BlenX language as the binding of the 
substrate and the enzyme is described as complex formation through specific binding sites of the boxes 
representing the proteins. In Figure [2} the types S and E are compatible and equipped with complexation 
and decomplexation rates. After complex-formation, S and E communicate and the internal behavior 
of the substrate box is changed into the behavior of the product (the ch(x,P) action modifies the type 
of the binder x into P). The new product has binding affinity no more to the enzyme and an infinite 
rate decomplexation occurs to release the enzyme E. Decomposition of the nonlinear term into ele- 
mentary reactions calls for the rate constants of each step. The rate constant of the catalytic reaction 
h = Vmax/[Etot] is easily obtained from the Michaelis-Menten formula. The dissociation constant of the 
ES complex ^2 = • — ^3 is obtained from the Michealis-Menten constants and it is supposed to be 
low because the ES complex is assumed to be stable. The rate constants of the reversible complexation 
(^1 and k2) can be chosen among several combinations by ensuring that the association rate is larger than 
the dissociation rate of the ES complex. Furthermore, we know that the catalytic step is the rate limiting, 
thus ki is chosen to be much larger than ^3. Our option determines the time of the simulation, thus values 
of the rate constants have to be carefully selected. In isolated systems we can scale down the constants 
easily in order to speed up the simulation. However the rates of the reversible complex formation have 
to be fast enough and cannot be limiting in a large model. 

The proper rate constants describing our compound Michaelis-Menten module have been selected by 
taking the minimum amount of substrate during the reaction and setting the initial (total) concentration of 
enzyme to Smin -0.1. As a consequence, we get ^1 = ^3 • 1000 = Vmax/[Etot] ' 1000 and = Km-ki—k^. 
Note that the selection of feasible parameters must lead to a positive value of ^2; and the total concentra- 
tion of enzyme has to be globally lower than the substrate with a large extent. 
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2.2 Simulation results 

We analyzed two models with a complex rate function and with an exact solution of the Michaelis- 
Menten kinetics. First, we converted the concentrations of the deterministic system into molecule num- 
bers through a transformation on the parameters using a scalar constant a defined as 1/ {Na • 10~^ • V), 
where Na is the Avogadro number and V is the volume of the modeled system (see the method of | 36 |). 
Then we set the models to different initial conditions for the substrate and run 200 simulations for each 
initials. The rates of product formation have been derived from the simulation results and these values 
are plotted over the initial amount of substrate molecules. It provides a saturation curve of the Michaelis- 
Menten kinetics (Figure [3]). 



v^g^= 60 #/min 










-—deterministic complex rate function (QSSA) 




—deterministic exact solution 




o stochastic complex rate function 




□ stochastic exact solution 



1000 2000 3000 4000 5000 6000 

Number of initial substrate molecules [#] 

Figure 3: Simulation results of the BlenX model fits the deterministic saturation curve well. At low 
amount of initial substrate molecules the original assumptions of Michaelis and Menten do not match 
and the exact solution diverges from the complex function. 

Calculation of the parameters is based on the assumptions shown previously. Km is set to 300# (# refers to 
number of molecules) and v^ax is 60#/mm. The total enzyme amount (\Etot |) is set to 60# molecules and 
the Michaelis-Menten constants define and the ratio of ki to ^2- The chosen parameters also satisfy 
the assumption that the value of ki (drrr'/ {min • mol)) is much larger than the value of ^2 (l/min). This 
condition may be suited by different rates of ki and k2, although these options only influence the speed 
of the reaction (and our simulation), but does not change the result (data not shown). 

Table 3: Parameters for the Michaelis-Menten module, a is set to 0.00167 during the simulations. 



Parameter names 


Parameter values 


Parameter units 


h 


200 -a 


1 / {min ■ #) 




99 


1 / min 


h 


1 


1 / min 


Km 


0.5 a 


# 


^max 


0.1 a 


#/min 


E 


0.1 a 


# 



We compared the deterministic and stochastic simulation results with the "unpacked" and complex mod- 
ules with a parameter set shown in Table |3] The module built up with complex reaction and the one with 
elementary reactions shows us a good accordance with each other and also with the deterministic scheme 
(Figure [3]). Simulation results of the BlenX model fits the deterministic saturation curve well, although 
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when the original assumptions of Michaelis and Menten do not match, the exact solution diverges from 
the complex function as it has been shown previously |[37]| . When the enzyme is in excess of the sub- 
strate, the solution of the unpacked model differ greatly from the packed version as the assumption made 
for the QSSA is not more valid for the system (Figure |3] and Figure |4]). This is one limitation of the 
compound function (also in deterministic simulations). 





10 20 40 50 10 20 30 40 50 10 20 30 40 SO 



Time (min) 

Figure 4: Simulation for different number of initial substrates: la-c: S = 30#; 2a-c: S = 60#; 3a-c: 
S = 599#. Average number of product is plotted for each time step. l-3a: stochastic unpacked versions 
shown with standard deviation from the mean. l-3b: complex rate functions with Kfn = 300#, v^ax = 
60#/mm where standard deviation is plotted. l-3c: comparison of unpacked (dashed) and complex 
(solid) stochastic simulation results made with BlenX. 

Arkin and Rao assumed 03 that the reactions are isolated and the amount of enzyme is fixed - but in 
complex network this assumption seems to be "weak". Enzyme concentration has to be much less than 
the substrate concentration and in e.g. oscillatory systems the substrate is changed over time. In those 
cases, the minimum value of the substrate has to provide the base of the calculation (see in Section 4). 
We emphasize that decomposition of Michaelis-Menten kinetics is not always necessary, but in a com- 
positional modeling framework it has to be available (as a library). Assumptions have to be checked and 
the decomposition might be especially useful for further extensions of the model (when hidden details 
are becoming important). For instance, when an inhibitor of the enzyme is present or two substrates of 
the same enzyme are introduced, details of the complex reactions have to be elucidated. In Section 2, we 
provided a brief description of a template for enzyme kinetics in BlenX with a parameter search based on 
basic mathematical calculus. Implementation of the template-library into the CoSBi Lab platform BTIl 
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might automatize the method of parameter estimation as it is a software including inference tools (KInfer 
ll421 ). Research proceed in that direction. 

3 Hill kinetics 

Sigmoidal response curves are usually described by a Hill function. Multiple reactions and several 
components that form complexes are hidden inside the Hill equation without making the details available. 
The average of a nonlinear function (e.g. Hill function) is generally found to differ from the function of 
the average (43], thus the proper way of handling these mathematical functions in biological models is 
crucial. The Hill equation assumes that n molecules of an entity (e.g. ligand) bind to a scaffold (e.g. 
receptor) simultaneously ||28]| and intermediate states do not occur. This is physically possible only if the 
number of ligands is equal to 1 (n = 1), but in most cases it is far from reality. 

Simultaneous binding of X molecules to Y can be described as F + nX ^ ^ YXn where ^/ is the rate 

h 

constant of the forward reaction and is of the backward. At equilibrium, the ratio of bound to total 
receptors is given by the Hill equation 

^ _ Bound _ YXn _ 
^'^^ ~ Total ~ Y + YXn ~ 

where the dissociation constant is J = k^/kf and n provides the number of ligands. The steepness of the 
transition of the sigmoidal curve depends on the number of ligands (n) and / provides the number of 
ligands at which half of the receptors (Y) are bound (Figure [5]). 




50 500 5000 

X 

Figure 5: Hill function for different Hill coefficients (n = 2, 3,4). When FHiii{X) is 0.5, / equals to 599#. 

Gene expression is known to be a particularly complex and noisy task ||7l[T0l[l2l[T3l[50|. Transcriptional 
regulation is often characterized by a sigmoidal Hill function. Nonlinearity arises from the assumption 
that the transcription factor forms multimers before binding to DNA (shown in Figure [6]), creating an 
abrupt switch between two states. The detailed mechanism behind the observed behavior is still unclear, 
but there are several hypothesis and models available |[29l[5Tl|52|. Difficulties of choosing a model has 
been proposed by ll53l . In the sequel, we investigate one simplified model of positive cooperativity that 
captures the requirements of containing only elementary reactions but still maintaining the sigmoidal 
property of the module. 
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3.1 Decomposition of the module 



Transcriptional factors (abbreviated as TFs) often form multimers during transcription ETi creating sig- 
moidal response of the system to the change of the transcriptional factors. Cooperativity widely occurs 
and the Hill equation is a good approximation of the underlying mechanisms, although it assumes simul- 
taneous binding of the TFs to the promoter region that is far from the realistic picture. As intermediate 
states have to occur during the reaction, sequential binding of the transcription factor to the promoter has 
been considered for this study. The following scheme approximates the Hill function when the interme- 
diate state (TF2) does not accumulate and positive cooperativity is present. Several other interactions 
are plausible ll29ll , but are not investigated in this article for the sake of simplicity. 



tf:TF tf:TF 
TF 



join , split 
-► 



join , split 



TF 



tf2:TF2 
_L 



TF2 



TF 



TF 



gtf2:GTF2 
_L 



GTF2_ 



TF2 




m:M 
J 



M 



5P/ 1 /\/\y 



Figure 6: Transcriptional regulation: Sequential binding of transcriptional factors {TFs) occurs then a 
homodimer (TF2) sits onto the promoter of a gene (G) and transcription results in messenger RNAs 
(mRNAs). Reactions are described through mass action kinetics with rate constants ki — ^4 and kms- 



BlenX offers a formal and efficient definition of join and split events of boxes as complex formation 
and decomplexation occur in biological systems. Transcription factors (TFs) form multimers (rF2, 3, 4) 
through a join event, enhancing the affinity of its binding onto the gene's promoter region (G). Positive 
cooperativity results in a low dissociation constant of the GTF2 trimer. The joined complex is able to 
transcribe mRNAs (Ms) of the gene through a split event that creates M boxes with the release of the 
active trimer. Bindings are assumed to be reversible. Elementary reactions of the system are summarized 
in Table II 

Sigmoidal curves are often measured in experiments providing specific Hill constants (such as n, J) 
to the Hill function. The Hill curve describing e.g. a transcriptional regulation scheme is not the proper 
way to apply stochastic calculations. Elementary steps of the sequential binding scheme contains four 
rate constants (^1, k2, ^3, ^4) that have to be determined from the constants n and /. Derivation of the 
missing parameters are calculated from the mass action kinetics describing the system. At equilibrium 
the intermediate complexes are assumed to be stable, thus 

TF2 ki , GTF2 h 

T = 1- and = — . 

TF^ k2 GTF2 k^ 

As the total amount of gene promoters do not change, G = Gtot — GTF2 leads to the term 

GTF2 _ TF^ 
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that is identical to a Hill function 

Bound _ X" 
Total 

Note that n = 2 and = (^2 • ^4)/(^i • ^3). The rate of transcription of M is • TF^/{J^ + TF^) = 
• {GTF2){Gtot) ^ where G^^r is a constant equals to 1 in this example. 

The Hill coefficient {n) and dissociation constant (/) determine only the relation of the four rate 
constants (/^ = (^2 •k/^)/{ki • k^)), but different values may satisfy the reaction scheme. We analyze 
several choices of ki,k2,k3, ^4 and describe them by the response coefficient (R) of the sigmoidal curve 
||55]| . The response coefficient allows us to measure the steepness of the transition in this reaction as it has 
been shown for other cases (such as Goldbeter-Koshland zero-order ultrasensitivity |55 |). R coefficient 
is defined as 5o.9/5o.i, the ratio of the signal (substrate) amount required to give 90% saturation relative 
to the amount required to give 10% saturation lISTTl . The dissociation constant (/) is chosen to be equal 
to 599# and the Hill coefficient (n) is 2. The complex reaction rate is TF^/ (599^ + TF^). In this case 
R = 9. 



Table 4: Reactions for Hill kinetics for the requirement of at least 2 ligands (TFs) 



Description 


Reactions 


Rate constants 


Units 




Homodimerization of TF 


2TF TF2 


ki 


1 / {min 


#) 


Dissociation of TF2 


TF2 2TF2 


ki 


l/min 




Formation of an active complex 


TF2 + G^GTF2 


k3 


1 / {min 


#) 


Dissociation of the active complex 


GTF2 ^ TF2 + G 


k4. 


l/min 




Synthesis of the mRNA 


GTF2 ^ M + GTF2 




l/min 





3.2 Simulation results 

We chose eight different set of parameters (see Table [s]) satisfying the following relation: = (^2 • 
k4){ki -ks). We measured the time average of the bound form (GTF2) in case of several levels of 
initial TF and calculated the response coefficient (R) and the actual Hill coefficient (n^) for each set of 
parameters. The derived Hill coefficient equals to log?>l/logR ll55]| . The derived dissociation constant 
(/) is calculated from the points as well as the root mean square error of the fit to the simulation results 
and the simulation point's error to the theoretical Hill function curve. 

If we compare the results of the complex function and the unpacked module, we see that when the 
assumption of Ki K2 is valid, the decomposed module gives a good fit to the theoretical Hill curve 
(Figure [7]). Our results agree with the observation of ll29l that for simple sequential binding schemes 
the only condition under which the Hill coefficient does accurately estimate the number of binding sites 
is when marked positive cooperativity is present. Furthermore, our analysis indicate that the larger the 
difference between the dissociation constants Ki and K2, the better the fit (e.g. compare set 1 to set 2). 

We also see that the Hill coefficient is not equal to the number of binding sites on the gene, but 
provides only a minimum value. 

The simulation results also shows behavior coincident to frequency modulation theory |[53l|56l where 
we see dense bursts of active transcription factors (GTF2) that results burstlike transcription giving 
similar downstream results that occur in concentration dependent transcription in deterministic models 
(data not shown). 
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Table 5: Multiple simulation results on the module of the Hill kinetics. Set refers to the deterministic 
version of the complex Hill function. Set 1-set 8 are different sets of parameters for the "unpacked" 
module. 
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TF 



Figure 7: The best fits (set 1,2,3,4) satisfy the assumption of positive cooperativity as Ki » K2 and 
provides a good fit to the theoretical complex function. 

4 A biological example 

After laying down the two widely used biological regulatory modules in BlenX, we move further and 
show how to use them during the composition of a biological case study. We chose circadian clock HTll 
as an example to present the arising noise of an oscillatory network containing complex rate functions. 
Circadian clock gives cells a daily rhythm to properly achieve several functions in a 24h periodic manner. 
This 24h oscillatory system is based upon a negative feedback loop producing a time delayed downregu- 
lation of a transcriptional activator. The presented scheme of the clock is adopted from a previous work 
done by f57l|. The model is a simplified model where nonlinearity has an important role in producing 
robust cycles, thus decomposition of complex rate functions is an interesting task. It is also an oscillating 
system where the assumptions have to stay valid for the whole model and where the search of param- 
eters is challenging. The proteins present in this system have several roles ll58]| . thus the possibility of 
choosing between complex rate functions or elementary reactions containing the required elements (e.g. 
an enzyme) of the reaction explicitly provides a more flexibile use of the language for later extension of 
the model. Compositionality remains an important and helpful advantage of BlenX in a template based 
environment. 

The system is divided into the following modules: (1) transcriptional regulation following Hill ki- 
netics (2) translation mechanism (assumed to follow mass action kinetics in this study) (3) homodimer- 
ization of clock proteins (CP) (4) formation of an inactive complex providing a negative effect inside 
the loop (5) There are three degradation term catalytically activated by enzymes (following Michaelis- 
Menten assumption) and the system also contains linear (so called background) degradation of the ele- 
ments. This network (Figure[8]) is built with complex rate functions and provides a 24h periodic oscillator 

the transcriptional step and the three degradation terms following Michaelis-Menten 
kinetics. We chose the number of enzymes having role in the system to be less than the corresponding 
substrates, making the assumptions of Michaelis-Menten kinetics valid. The parameters originating 
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Figure 8: Circadian clock model composed of transcriptional and enzyme catalyzed degradation mod- 
ules. 



from the complex functions are also scaled up to be fast enough. Thus the reactions assumed to be in 
equilibrium do not limit the system. The products of the enzymatic reactions are degraded immediately 
(with an infinite rate) after their production in order to serve the catalyzed degradation scheme in the 
system. 

We simply merge the modules and insert the boxes (enzymes and intermediates) of the novel entities. 
We also replace the events corresponding to the complex functions for the ones from the "unpacked" 
modules. This method can be easily automatized as it does not require the modification of the reac- 
tions that are independent of the substituted complex functions and the functions calculating the rate 
of complex reactions does not involve binders. The novel internal behavior of the boxes can be easily 
parallelized with the original ones. The composition of modules is shown in Figure [9] Simulation of 
the "unpacked" system (Figure [T0](1 -2a)) shows larger noise than the original (Figure [T0](1 -2b)), but 
still produces regular oscillations. It has been shown in several works |[32l|59l that with process calculus 
based languages dynamic models can be constructed and existing continuous models can be transferred 
into the stochastic framework providing additional predictions of the biological system results to the ex- 
isting models |38|. Herein we have to remind the reader that generally distributed reaction times have 
been also implemented into the BlenX framework recently [ 60]. It provides choices of the reaction time 
distribution for the stochastic simulation algorithm of Gillespie. In this way, abstracted rate laws can 
be handled stochastically that leads to a better quantitative tool for matching wet-lab experiments and 
in-silico results without breaking down the complex reactions into elementary steps. The use of this ex- 
tension fits well the idea of a template based modeling framework as, depending on the question the user 
asked, biological models might be characterized through complex rate laws and handled by generalized 
distributions of time; while templates (including only elementary steps) offer a straightforward, flexible 
and more precise way of compositional modeling in BlenX. 



5 Conclusion 

Two basic complex rate functions that are widely used in biological modeling have been decomposed 
into elementary reactions in the BlenX framework. The work proposed in this article provides a proper 
structure for Gillespie's algorithm. Furthermore, the templates offer a straightforward method for compo- 
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let 


CP : 


bproc 




#(cp,CP) 


[ nil ]i 


let 


CP2 : 


bproc 




#(cp2jCP2) 


[ nil ]j 


let 


IC : 


bproc 




#(ic,IC) 


[ nil ]j 


let 


TF : 


bproc 




#(tf ,TF) 


[ nil ]; 


let 


M : 


bproc 




#(mjM) 


[ nil ]; 


let 


CP : 


bproc 


= 


#(cp,CP) 


[cp!().ch(cp,CPUl_D).nil]; 


let 


CP2 : 


bproc 


= 


#(cp2,CP2) 


[cp2!().ch(cp2jCP2U2_D).nil]; 


let 


IC : 


bproc 


= 


#(ic,IC) 


[icl().ch(ic,ICU3_D).nil]; 


let 


TF : 


bproc 


= 


#(tf,TF) 


[ nil ]; 


let 


M : 


bproc 


= 


#(mjM) 


[ nil ]; 


let 


Ul : 


bproc 




#(uljUl) 


[ rep ul?().nil ]; 


let 


U2 : 


bproc 




#(u2jU2) 


[ rep u2?() .nil }; 


let 


U3 : 


bproc 




#(u3jU3) 


[ rep u3?() .nil ]; 


let 


CPU1_D : 


bproc 




#(cpjCPUl_D) 


[ nil ]i 


let 


CP2U2_D 


bproc 




#(cp2,CP2U2_D) 


[ nil ]; 


let 


ICU3_D : 


bproc 




#(iCjICU3 D) 


[ nil ]; 


let 


G : 


bproc 




#(g>G) 


[ nil ]; 


let 


TF2 : 


bproc 




#(t2f,TF2) 


[ nil ]; 


let 


GTF2 : 


bproc 




#(gtf2jGTF2) 


[ nil ]; 



(b) 



(c) 



synthesis ofM 

when(M: iTRANSCRIPTION) new(l); 

where 

let TRANSCRIPTION : function = 



(kms*pow( |TF I jn))/(pow(3jn) + pow( | TF | , n) ) ; 



synthesis ofM 

when(TFjTF: :rate(kl)) join<TF2); 
when(TF2: :rate(k2)) split(TF^TF); 
when(TF2,G: :rate(k3)) join(GTF2); 
when(GTF2: :rate(k4)) split (GjTF2)j 
when (GTF2 :: TRANSCRIPTION) split (GTF2,M); 
where 

let TRANSCRIPTION : function = kms*|GTF2|; 



(d) 



degradations of CP, CP2 and IC through Michaelis-Menlen kinetics 
when(CP: :CP_DEG) delete(l)j 
when(CP2: :CP2_DEG) clelete(l)j 
when(IC: :IC_DEG) split(NiljTF); 

where 

let CP_DEG : function = kpl* | CP | /(Kml+ | CP | ); 
let CP2_DEG : function = kp2* | CP2 | /(Km2+ | CP2 | ) j 
let IC_DEG : function = kp3* | IC | /(Km3+ | IC | ); 



degradations of CP, CP2 and IC through Michaelis-Menlen kinetics 
when(CPUl_D: :inf) delete(l)i 
when(CP2U2_D: :inf) delete(l); 
when(ICU3_D: :inf) split(NiljTF); 

(UljCPjrate(kpll)jrate(kpllr)jrate(kpl2))j 
(UljCPUl_Dj0>inf,0), 

(U2jCP2,rate(kp21),rate(kp21r),rate(kp22)), 
(U2 J CP2U2_D,0j inf J 0) j 

(U3jICjrate(kp31)jrate(kp31r)jrate(kp32))j 
(U3jICU3_Dj0,inf,0) 



(f) 



degradation ofM 

when(M: :rate(kmd)) delete(l)j 
translation of M into CP 
when(M: :rate(km2)) split(MjCP)j 
homodimerization of CP and dissociation of CP2 
when(CPjCP: : rate(ka)) join(CP2); 
when(CP2: :rate(kd)) split(CPjCP); 
inactive complex formation and dissociation 
when(CP2,TF: :rate(kica)) join(IC)j 
when(IC: :rate(kicd)) split(CP2,TF); 
degradation terms of CP, CP2 and IC 
when(CP: :rate(kcpd)) delete(l)j 
when(CP2: : rate(kcp2d) ) delete(l); 
when(IC: :rate(kcp3d)) split(Nil,TF); 



(g) 



Figure 9: BlenX source code. The boxes of the original model are shown in (a) while the "unpacked" 
version is in (b). Composition of (a) and (b) is a straightforward job by parallelization. Events of the 
original model are in (c), (e), (g), while (d), (f), (g) contains the "unpacked" version of the model. Note 
that there is no change in (g). The substituited modules are highlighted (bold font) in the text. Parameters 
are provided upon request. 
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Figure 10: Simulation results for the stochastic model containing complex rate functions (l-2b) and 
for the "unpacked" versions (l-2a) in case of a conversion factor a=0.000167 (la-b) and a=0.0000167 
(2a-b). The total amount of CP in the system is plotted as solid curves; dashed curves represent the 
messenger (M) while dotted points demonstrate the free transcription factors (TF) in the model. 



sitional modeling. The usage of templates provides a less error prone method in modeling as assumptions 
are taken into account during the model composition. We exemplified our approach on a circadian clock 
model. On the top of the first set of predefined templates, a higher level of abstraction and higher level 
of model composition might be exploited in the future by building larger templates. For instance, a bi- 
ological switch might be built from two Michaelis-Menten functions, often called Goldbeter-Koshland 
module |[55l . Although the simulation time increases by decomposing complex reactions into elemen- 
tary ones, we think that the modeling advantage is worth pursuing. Optimizations and simplification can 
be postponed to the model of the overall system. We will continue investigating the decomposition of 
complex reactions in order to build a library that will allow drag-and-drop style of modeling complex 
behavior relying on predefined bricks. 
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